

#libraries and modules
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.interpolate import UnivariateSpline
from scipy.stats import gaussian_kde, beta, lognorm
from lmfit import minimize, Parameters,report_errors

def betadist(betaparams,B,pcF):
    #defining beta distribution parameters
    a=float(betaparams['a'].value)
    b=float(betaparams['b'].value)
    loc=float(betaparams['loc'].value)
    scale=float(betaparams['scale'].value)
    #creating fitted data
    model_pcF=beta.pdf(B,a,b,loc=loc,scale=scale)         
    #returning residual
    return (model_pcF-pcF)
    
def betadist_leastsquare_fitting(B,Bpdf,ai,bi,loc,scale):
    #defining power law variables
    betaparams = Parameters()
    betaparams.add('a', value=ai, vary=True, min=0.001, max=None)
    betaparams.add('b', value=bi, vary=True, min=0.001, max=None)
    betaparams.add('loc', value=0, vary=True, min=-5, max=0)
    betaparams.add('scale', value=90, vary=True, min=70, max=120)
    #minimizing residuals
    result = minimize(betadist, betaparams, args=(B,Bpdf))
    Bfit = np.linspace(betaparams['loc'].value,betaparams['scale']+betaparams['loc'].value,90*10+1)
    Bpdffit=beta.pdf(Bfit,
                    betaparams['a'].value,
                    betaparams['b'].value,
                    loc=betaparams['loc'].value,
                    scale=betaparams['scale'].value)
    
    return Bfit,Bpdffit, np.sqrt(np.mean((result.residual)**2)),betaparams['a'].value,betaparams['b'].value,betaparams['loc'].value,betaparams['scale'].value




def grad_dist(curdata,ax):
    #gaussian kernel density estimation           
    x=np.linspace(0,90,91)
    kde=gaussian_kde(curdata)
    line,=ax.plot(x,kde(x), '-', label=hl[j][0:-2]+"0")
    curcolor=plt.getp(line,'color')

##    #plotting histogram
##    ax.hist(curdata,bins=int(max(curdata))+1,
##             normed=True, histtype='step',
##             color=curcolor, linewidth=0.5)

    #defining splines for kd function and corresponding 1st and seconde derivatives
    x=np.linspace(0,90,901)
    s=UnivariateSpline(x,kde(x),s=0,k=3)
    s1=UnivariateSpline(x,s(x,1),s=0,k=3)
    s2=UnivariateSpline(x,s1(x,1),s=0,k=3)

    #identifying local maxima (where s1=0, s2<0, and s>0.005)
    maxima=s1.roots()[np.where(s2(s1.roots())<0)[0]]
    maxima=maxima[np.where(s(maxima)>0.005)[0]]
    #s_max=maxima[-1]
    s_max=maxima[np.argmax(s(maxima))]
    ax.plot(s_max, s(s_max),'o', color=curcolor)

    #identifying steepest segment after maxima (where x>=maxima, s1<0)
    x2=x[np.where(x>=s_max)[0]]
    slope=s1(x2)
    s1_min=x2[np.argmin(slope)]
    ax.plot(s1_min,s(s1_min),'o', color=curcolor)

    print round(s_max,1),round(s1_min,1)

    return round(s_max,1),round(s1_min,1)

def power(powerparams,hl,B):
    hl=np.array(hl)
    B=np.array(B)
    #defining power distribution parameters
    a=float(powerparams['a'].value)
    b=float(powerparams['b'].value)
    #creating fitted data
    model_B=a*(hl**b)
    #returning residual
    return np.array((model_B-B))
    
def power_leastsquare_fitting(hl,B,ai,bi):
    #defining power law variables
    powerparams = Parameters()
    powerparams.add('a', value=ai, vary=True, min=0.001, max=None)
    powerparams.add('b', value=bi, vary=True, min=None, max=-0.00001)
    #minimizing residuals
    result = minimize(power, powerparams, args=(hl,B))
    hlfit = np.linspace(30,110,100)
    Bfit=powerparams['a'].value*(hlfit**powerparams['b'].value)
    
    return hlfit,Bfit, np.sqrt(np.mean((result.residual)**2))






#MAIN
plt.close('all')
np.set_printoptions(suppress=True)




lith=['CC','FC','LS','UM']
lithnum=[2,3,4,5]
hl=['n3gr','n5gr','n7gr','n9gr','n11gr']



grad=pd.read_csv('lith_fa_gr_3_5_7_9_11.csv',sep=',',names=['x','y','lith','fa','n3gr','n5gr','n7gr','n9gr','n11gr'],index_col=[1,0])
grad=grad.reset_index(drop=True)
grad.n3gr=np.round(grad.n3gr,1)
grad.n5gr=np.round(grad.n5gr,1)
grad.n7gr=np.round(grad.n7gr,1)
grad.n9gr=np.round(grad.n9gr,1)
grad.n11gr=np.round(grad.n11gr,1)
grad.fa=np.round(grad.fa,4)
grad=grad.dropna(axis=0, how='any')


fig,ax=plt.subplots(nrows=len(lith), ncols=1, sharex=True, sharey=True)
df2=pd.DataFrame(data=None,index=None, columns=['Li','hl','fmax','f1min','25','50','75','100'])
    
for iL in range(len(lith)):
    
    
    plt.sca(ax[iL])
    if iL==0:
        curgrad=grad[(grad.lith<=lithnum[iL])]
    else:
        curgrad=grad[(grad.lith==lithnum[iL])]

    #sorting dataframe according to 'fa'
    curgrad=curgrad.sort('fa')
    #removing "ridge" pixels (fa<5 pixels)
    #curgrad=curgrad[(curgrad['fa']>=5)]
    #resetting index
    curgrad=curgrad.reset_index(drop=True)

   
    
    for j in range(len(hl)):
        #computing expanding mean for gradient at current scale
        curgrad['cum_mean']=pd.expanding_mean(curgrad[hl[j]])
        i=curgrad['cum_mean'].idxmax(axis=0)

        #computing expanding mean for gradient at current scale
        curgrad['cum_max']=pd.expanding_max(curgrad[hl[j]])
        imax=curgrad['cum_max'].idxmax(axis=0)

        if imax>i:i=imax
        maxfa=curgrad['fa'].values[i]
        minfa=curgrad['fa'].values[i]/2.

        if minfa<5:minfa=5
        if maxfa<10:maxfa=10
##        print "FAmaxgrad: ", curgrad['fa'].values[i]
##        print "minFA: ", minfa
        
        #removing pixels with 'fa'> 'fa' at first occurrence of maximum gradient
        xcurgrad=curgrad[(curgrad.fa<=maxfa) ]
        xcurgrad=xcurgrad[(xcurgrad.fa>=5) ]
        f_max,f1_min=grad_dist(xcurgrad[hl[j]].values,ax[iL])
        
        dfrow={'Li':lith[iL],'hl':int(hl[j][1:-2]+'0'),'fmax':f_max, 'f1min':f1_min}

        df2=df2.append(dfrow,ignore_index=True)
        curax=ax[iL]
        curax.set_ylabel(lith[iL])
    
    curax=ax[-1]
    curax.set_xlabel("gradient, degrees")
    plt.legend(fontsize='x-small')

print df2
df2.to_csv('gradient_inventory_analysis_results.csv',sep=',')






fig,ax=plt.subplots(nrows=2,ncols=1, sharex=True, sharey=True, figsize=(3.5,6))
main_plot=fig.add_subplot(111, frameon=False)
main_plot.set_ylabel('slope gradient, degrees\n\n ')
main_plot.set_xlabel('\n\nbase length, m')
main_plot.set_xticklabels([])
main_plot.set_yticklabels([])
main_plot.set_xticks([])
main_plot.set_yticks([])
for j in range(len(lith)):
   
    
    cur_df2=df2[(df2.Li==lith[j])]
    hl=cur_df2.hl.astype(float).values
    f=cur_df2.fmax.astype(float).values
    f1=cur_df2.f1min.astype(float).values

    plt.sca(ax[0])

    ai,bi=30,-1

    line,=plt.plot(hl,f, '.')
    curcolor=plt.getp(line, 'color')
    
    
##    hl2=np.linspace(hl.min(),hl.max(),100)
##    s1=UnivariateSpline(hl,dydx,s=1,k=2)
##    dydx=s1(hl2)
##    plt.plot(hl2,dydx, ':', color=curcolor)

    hl2,s,RMSE=power_leastsquare_fitting(hl,f,ai,bi)
    plt.plot(hl2, s, '-', color=curcolor, label=lith[j])
    
    if j==len(lith)-1:
        plt.title('maximum pcF gradient', fontsize='small')
        plt.legend(fontsize='x-small', loc='upper right')
        ax[0].set_yticks((15,30,45,60,75))

    plt.sca(ax[1])
    plt.plot(hl,f1, '.', color=curcolor)

##    s2=UnivariateSpline(hl,d2ydx2, s=1,k=2)
##    d2ydx2=s2(hl2)
##    plt.plot(hl2,d2ydx2, ':', color=curcolor)

    hl2,s1,RMSE=power_leastsquare_fitting(hl,f1,ai,bi)
    plt.plot(hl2,s1, '-', color=curcolor)

    if j==len(lith)-1:
        plt.title('maximum pcF curvature', fontsize='small')
        ax[1].set_xticks((30,50,70,90,110))
        ax[1].set_yticks((15,30,45,60,75))
        ax[1].set_ylim(15,75)
        ax[1].set_xlim(25,115)

fig.tight_layout()









plt.show()

        


            
            

            



                




   


   
            


